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Abstract 



In this paper we study the Lyapunov stability and Hopf bifurcation in a 
biological system which models the biological control of parasites of orange 
plantations. 
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1 Introduction of the Mathematical Model 



In this work we study a system of four coupled differential equations ([T]) 
which models the interaction between two biological species, each presenting 
two stages in their metamorphosis, living in a common habitat with limited 
resources. 

The differential equations analyzed here are 

P' = ^ = 01 f 1 - — ^ M - (ai + /3i)P - hPG 



dt \ C\ 

M' = ^ = aiP-fiiM (1) 

= ^ = (1^2(1 --^G-{a2 + P2)L + k2PG 
dt \ C2J 

dG 

G' = — — = a2L — 112G. 
dt 

This model — an elaboration of Lotka-Volterra equations, taking into ac- 
count the stages or compartments in the biological populations — was pro- 
posed bv Yang and Ternes [U, Ej and Ternes [3] for a study of the biological 
controU of orange plantations leaf parasites P, which is a pre-adult stage for 
M, by their natural enemies L, which is an early stage for G. 

Other differential equations have been proposed as models for interacting 
populations partitioned in compartments, representing several situations of 
biological interest. See, among many others, Hethcote et al. [1], Jacquez and 
Simon [5] and Godfray and Waage [6]. 

In [U E] and [3] P and M are the densities of pupaj^ and female adults 
of Phyllocnistis citrella (which in its larval^ stage is the citrus leafmineJl), L 
and G are the densities of larvae and female adults of its native parasitoid 



^http:/ /en. wikipedia.org/wiki/BiologicaLcontrol 
^http: / / en.wikipedia.org/ wiki/Pupa 
'^littp://en. wikipedia.org/wiki/Larva 
^Lttp:/ /www.agrobyte.com.br/minadora.htm; http://en.wikipedia.org/wiki/Citrus 
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Galeopsomyia fausta (whose larvae feed on the pupae 

The meaning of the parameters in ([1]), where the notation of [U |3] has 
been preserved, is as follows: ai is the rate of pupae that give rise to adults 
M, j3i is the mortality rate of pupae, /ii is the mortality rate of adults M, 
01 is the rate of eggs that give rise to pupae, Ci is the carrying capacity of 
the population M, a2 is the rate of larvae that, evolving through pupae, give 
rise to adults G, [32 is the mortality rate of larvae and pupae, ^2 is the rate 
of mortality of adults G, 02 is the ovipositioij^ rate of the parasite and C2 is 
the carrying capacity of the population G. Here we assume that the pupa 
(respectively larva) population decreases (respectively increases ) at a rate 
proportional to P — G encounters that is kiPG (respectively k2PG). 

This model represents the evolution of female populations. If necessary 
the male populations can be estimated using the sexual ratio of each species. 

Remark 1.1 All the parameters ai, /ii, 0i, Ci, /ci, 02, /?2, At2, 02, C2, ^2 are 
positive. As the damage to the P population must be larger than the benefit 
to the L population it is natural to assume that ki> k2. 

Here will be established the location and the stability character of the 
equilibria of ([T]), four in number. Also is determined the bifurcation variety 
in the space of parameters, representing the transition from asymptotically 
stable to saddle type at the equilibrium point with positive coordinates, 
representing the coexistence of the two species. See Theorem 12.41 and its 
Corollary [231 

Fixing the all the parameters in ([1]) to biologically feasible values, taken 
from [3] and [7], but letting the interaction coefficients ki and k2 vary in a 
positive quadrant, the nature of the bifurcation phenomenon in this plane 
by crossing the bifurcation curve is established. See Theorem 13.11 and Figure 
[U This is done by means of a computer assisted calculation of the first 

'^http://www. seea.es/conlupa/AlbertoWeb/framesparasitoides. htm. This site has im- 
pressive photos of hosts and parasitoids. 
^en. wikipedia.org/wiki/Oviposition 
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Lyapunov coefficient, found to be positive. Tlie Hopf bifurcation analysis 
of this point implies that the bifurcating periodic orbit is asymptotically 
unstable, of saddle type which surrounds an attracting equilibrium with small 
attracting basin. The dependence of the bifurcation curve on the parameter 
C2 is studied in Theorem 13.31 and illustrated in Figure O 

In Section m the implications of the results in this paper are discussed and 
interpreted from a wider perspective. 

2 Stability Analysis of Equilibria 

Assume the following notation: 

The differential equations ([1]) have four equilibria 

.4i = (Pi,Mi,ii,Gi) = (0,0,0,0), (3) 

A, = = (l - ^) (l - i) .m) . (4) 

^3 = (P. Ma. L„ G,) = (o. 0. 5J (l - i-) . (l - i-) ) . (5) 

and 

A= (P4,M4,L4,G'4), (6) 

where 

= ^ ---^-r- I ^^0^ / ^ _ 1 j _ ^^c2/i;i (1 - 1^ 





102 + filciC2kik2 




ciai02 




'102 + ^l\ciC2kik2 




C2/i2ai01 




0102 + filciC2kik 




02^101 



-^4 = ^^2 J. .r'"T.".2^ ^ \ ( ClyUlfc2 ( 1 - + «102 (^^ ~ ^ 

G4 = 2 J, ^ — — - ( ci/ii/i;2 ( 1 - 4- I + «i02 (1 - ^1 

ai0102 + /iiCiC2fclfc2 V V RlJ V -R2/, 
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Remark 2.1 If Ri > 1 and R2 > 1, then the equilibria Ai, A2 and A^, have 
only non-negative coordinates. If ki < where 

Q;i0i (1 - ^) 

klmax ~ 7 C^) 

C2/il [^-R-,) 

then the coordinates of the equilibrium ^4 are also non-negative. 
The Jacobian matrix of ([T]) at x = (P, M, L, G) G has the form 



J(x) 



a\ — /ii 

k2G -a2-/32 02-^ + A;2P 

V a2 -Ai2 / 



(8) 



while its characteristic polynomial is given by 

p(A) = c/et(J(x) - A/) = 61 02 + a2(/ii + X)kik2PG, (9) 

where 

Gi = (/ii + A)(ai + /3i + fciG + A) - ai (01 - ^^i^ 

V Cl 

and 

82 = (/i2 + A)(a2 + (32 + X)-a2 (02 - + A;2P ) . 



C2 

Recall that an equilibrium point xg is said to be a saddle of type n — p ii 
the Jacobian matrix J(xo) has n eigenvalues with negative real parts and p 
eigenvalues with positive real parts. 

Theorem 2.2 If Ri > 1, P2 > 1 and ki < ki^^^ then: 

1. The equilibrium Ai is a saddle of type 2-2; 

2. The equilibrium A2 is a saddle of type 3-1; 

3. The equilibrium A^ is a saddle of type 3-1. 
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Proof. From the eigenvalues of J{Ai) are given by 



Ai = + A + /ii) + I + A + + 4ai0i[l - 



A2 = + A + /^i) - + A + /^i)' + 4ai0i[l - 



A3 = -|(a2 + /52 + /i2) + I J ("2 + /32 + /^2)2 + 4a202[l - 4 



A4 = -|(a2 + /52 + /i2) - |\/ ("2 + /32 + /i2)^ + 4a202[l - 7^], 



and satisfy: Ai>0, A2<0, A3>0 and A4 < 0. This proves the first 
assertion. 

From ([9]) the eigenvalues of J (^2) are given by 



Ai = -^(ai + A + /ii) + (ai + A + - 4«i0i[l - 



A2 = -^("1 + A + yui) - ^ Y ("1 + Pi + fj^iy - 4ai0i[i - 

A3 = -5(a2 + /52 + Ai2) + 



i^(«2 + /?2 + /i2)2 + 4a202[l - Tfe] + 4^[1 - ir]^2, 
A4 = -|(a2 + /52 + /i2)- 

yia, + + /i2)2 + 4a202[l - + 4^[1 - ^Jfc^. 
Is immediate to see that A3 > and A4 < 0. If 

4ai 

then Ai and A2 are complex with negative real parts and if 

1 



^1 < 7^[(tti + A + /ii)^ + 4/ii(ai + (3i)] 
then Ai < and A2 < 0. This proves the second assertion. 



From the eigenvalues of J {A3) are given by 

Ai= -liai + f3i + fii + C2h[l - ^])+ 

l^{ai + A - /xi + C2/ci[l - + 4ai0i, 



i + A - /ii + C2fci[l - + 4ai0i 



1 1 / 1 

A3 = («2 + /52 + /U2) + o\/ («2 + /?2 + /i2)^ - 4^202 [1 " ^] 

z z V ^2 



1 1 / 1 

A4 = --(^2 + ^2 + /i2) - 2 Y ("2 + /32 + 1^2)^ - 4a202[l " -^] 

It follows that Ai > 0, A2 < 0. If 

(t>2 > 7^[(«2 + /52 + /i2)^ + 4/i2(a2 + /?2)] 

4q;2 

then A3 and A4 are complex with negative real parts and if 
02 < T^[(a2 +P2+ /i2)^ + 4/i2(a2 + P2)] 

4^2 

then A3 < and A4 < 0. This proves the last assertion. 



For the sake of completeness we state the following lemma which is a 
particular case of the Theorem of Routh-Hurwitz. See [8], p. 62. 

Lemma 2.3 The polynomial L{X) = oqA^ + aiX^ + a2A^ + a3A + a^, ao > 0, 
with real coefficients has all roots with negative real parts if and only if the 
numbers 01,02,03,04 are positive and the inequality 

A = oi 02 03 — oq O3 — o^ 04 > 

is satisfied. 
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Theorem 2.4 // i?i > 1, i?2 > 1 and ki < ki^^^ then all the coefficients of 
the characteristic polynomial of JIAa) are positive. Therefore, if 

A = ai a2 as — 03 — al a4, > 0, (10) 

where 

ai = ai + Pi + fii + a2 + (32 + + ^1^4, 
0,2 = ^^i^M4 + ^^^6*4 + («! + Pi + ni + kiGi){a2 + P2 + 1^2) ■, 

Cl 02 

0-3 = ("1 + A + Ail + ^1(^4) ^4 + ("2 + /52 + yU2) M4 + a2kik2PAGi, 

C2 Cl 

04 = a2ai</.i (k2 f 1 - ) P4 + — ( 1 - 1^ I ^4 ) , 



\ \ Rl J Cl \ i?2 / / 

then the differential equations (QP have an asymptotically stable equilibrium 
point at A4. If 

A < 

then A4 is unstable. 

Proof. From the characteristic polynomial of J(^4) is given by 



[A^ + («i + + kiG^)\ + ^M4][A2 + («2 + /52 + ^^2)\ + ^6-4] 

Cl C2 
+a2/Ul^l ^2-^4^4 + a2^1^2-P4G4A, 

which can be written as 
A^ + \^[ai + (3i + + a2 + (32 + ^^2 + ^1^4] 

+A2[^ii^M4 + + («i + A + /ii + A;iG'4)(«2 + P2 + ^^2)] 

+A[(«i + A + /ii + A;iG4)^G4 + (a2 + /32 + 7^2) ^M4 + a2A;ifc2P4G'4] 

+ ^M4^G'4 + «2/xiA;iA;2P4G'4. 

Now it is simple to see that the coefficients of the characteristic polynomial 
are given by 01,02,03,04 above. From the hypotheses these coefficients are 
positive. The stability at follows from Lemma [2.31 
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The following corollary is immediate from the fact that aj > 0. 

Corollary 2.5 The Jacobian matrix J{A4) has a pair of complex eigenvalues 
with zero real part if and only if 

al — Oi 02 03 + al = 0, (11) 

where ai are defined in Theorem\2.4 



In next section we study the stability of under the condition (fTTj) . 
complementary to the range of validity of Theorem I2.4[ 

3 Hopf Bifurcation Analysis 
3.1 Generalities on Hopf Bifurcations 

The study outlined below is based on the approach found in the book of 
Kuznetsov [9], pp 175-178. 

Consider the differential equations 

x' = /(x,/i), (12) 

where x G and /i G is a vector of control parameters. Suppose (fT2!) 
has an equilibrium point x = xq at /i = yUo and represent 

F(x) = /(x,/xo) (13) 

as 

F(x) = Ax+ ^ fi(x,x) + \ C(x,x,x) + 0(||x||^), (14) 
where A = /x(0,/io) and 



a(x, y, z) = J2 J. .-n^^- 



,,,=1 di, dii 



for z = 1, 2, 3, 4. Here the variable x — xq is also denoted by x. 

Suppose (xo, /io) is an equilibrium point of (fT2|) where the Jacobian matrix 
A has a pair of purely imaginary eigenvalues A3_4 = ±ic<Jo, cuq > 0, and no 
other critical (i.e., on the imaginary axis) eigenvalues. 

Let p, g G C"^ be vectors such that 

4 

Aq = iuo q, A^p = -iujQ p, {p, q) = ^Piqi = 1. (17) 

i=l 

The two dimensional center manifold can be parameterized by w G = 
C, by means of x = H{w, w), which is written as 

H{w,w) = wq + wq + hjkW^w'' + 0{\w\^), 

2<j+fc<3 •' 

with hjk G C^, hjk = hkj. 

Substituting these expressions into f|T2l) and (|T4l) one has 



H^{w, w)w' + Hisiw, w)w = F{H{w, w)). (18) 

The complex vectors hij are to be determined so that equation f|T8|) writes 
as follows 

w' = iujQW + - G2iw\w\^ + 0(|w|^), 

with G21 G C. 

Solving the linear system obtained by expanding (JTSil . the coefficients of 
the quadratic terms of (fT3l) lead to 

h,, = -A-^B{q,q), (19) 

h2^ = {2iuJoh-A)-^B{q,q), (20) 
where I/^ is the unit 4x4 matrix. 
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The coefficients of tlie cubic terms are also uniquely calculated, except 
for the term w'^w, whose coefficient satisfies a singular system for h2i 

{tuoh - A)h2i = C{q, q, q) + B{q, /iso) + 25(g, /in) - Gsig, (21) 

which has a solution if and only if 

(p, C(g, g, q) + h^o) + 25(g, /m) - Gsig) = 0. 

Therefore 

G21 = {p,C{q,q,q)+B{q,{2zuJoh-A)-'B{q,q))~2B{q,A-'B{q,q))), (22) 

and the first Lyapunov coefficient li - which decides by the analysis of third 
order terms at the equilibrium its stability, if negative, or instability, if pos- 
itive - is defined by 

/i = -^ReG2i. (23) 

A H op f point (xo,/io) is an equilibrium point of (fT2|) where the Jacobian 
matrix A has a pair of purely imaginary eigenvalues A3^4 = iicuo, ujq > 
and no other critical eigenvalues. At a Hopf point, a two dimensional center 
manifold is well-defined, which is invariant under the fiow generated by (fT2l) 
and can be smoothly continued to nearby parameter values. 

A Hopf point is called transversal if the curves of complex eigenvalues 
cross the imaginary axis with non-zero derivative. 

In a neighborhood of a transversal Hopf point with /i 7^ the dynamic 
behavior of the system f[T^ . reduced to the family of parameter-dependent 
continuations of the center manifold, is orbitally topologically equivalent to 
the complex normal form 

w' = {'y + iuj)w + liw\w\'^ ^ (24) 

w G C, 7, and li are smooth continuations of 0, uoq and the first Lyapunov 
coefficient at the Hopf point 0, respectively. When /i < {li > 0) a family 
of stable (unstable) periodic orbits can be found on this family of center 
manifolds, shrinking to the equilibrium point at the Hopf point. 
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3.2 Hopf Bifurcation in the Biological Model 

In this subsection we analyze the stabihty at A4 given by under the 
condition (fTTj) . From (fT2|) write the Taylor expansion (HM of /(x). Thus 

f ~{ai + P,) - hG, 0i(l-^) ~hP^ \ 



A 



cti — /ii 

k2G^ -ia2+f32) 02(1-^) + ^2^4 

V ^2 -Ai2 / 

(25) 

and, with the notation in f[T^ to (ITB]) . one has 

F(x) - Ax = f-^li^ - kiPG, 0, + k2PG, 0] . (26) 

V Cl C2 / 



From (111, ([T5]), ([T6]) and ([26]) one has 

5(x,y) = (5i(x,y),0,53(x,y),0), (27) 



where 



5i(x, y) = -— X2 1/2 - /i;i(a;i 1/4 + yi), 

Cl 

53(x, y) = -— Xi 2/4 + k2{xi y4 + X4 yi), 

C2 



and 



C(x,y,z) = 0. (28) 
To pursue the analysis consider the following table of specific parameters 



ttl 


= 0.7 


Pi 


= 0.003 


/ii 


= 0.6 


01 = 2.3 


Cl 


= 400000 


a2 


= 0.3 


P2 


= 0.0015 


fJ'2 


= 0.4 


02 = 4 


C2 


= 100 



(29) 



taken from ^ and [3], where their biological feasibility in Brazilian fields is 
discussed. 

With the above parameter values the differential equations ([1]) are in fact 
a two parameter system of differential equations where the parameters are 
ki and k2 and can be written equivalently as 

x' = /(x,A;i,A;2), (30) 
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with /(x, ki, defined by the right-hand sides of ([T]). 

With the parameter values of table fl29|) . the equilibrium point A4 has 
the following coordinates 

800000 (1.425 - 64.764 ki) _ 933333.333 (1.425 - 64.764 ki) 



4.508 + 1.444- 10^ A;i A;2 ' 4.508 + 1.444 • 10^ fci fcs 

444.444 (1.216 + 85550.400 ^2) „ 333.333 (1.216 + 85550.400 A^a) 

(jrA 



4.508 + 1.444 ■ 10^ ki k2 4.508 + 1.444 ■ 10^ ki fca 

while -Ri, i?2 and given by ([2]) and ([7j), have the form 

Ri = 3.81697, R2 = 9.95025, A;i_ = 0.0220159. (31) 

From the above equation and the Remark 11.11 the set of admissible parame- 
ters is given by (see Fig [T]) 

S = {{ki, k2)\0<ki< = 0.0220159 and < A;2 < k^}. (32) 

In this set 5 the curve S = A~^(0) is well-defined (see (|TT1) ). where A is 
given by 

1699.422 - 233762.372A;i - 6.860 ■ 10^fc2 - 1.114 ■ lO^A;^ - 2.175 ■ lO^^fc^ - 
4.994 ■ 10^°A;iA;2 - 2.147 ■ lO^kf - 4.079 ■ lO^^k^ - 1.529 ■ lO^^kikj - 
7.809 • 10^^klk2 - 4.319 • lO^'^klkj - 1.540 ■ lO^^kik^ - 4.755 • 10^^A;fA;2 - 
1.752 ■ lO^^kfkl - 6.741 ■ 10^^A;^A;^ + 2.940 ■ lO^^k^k^ - 1.703 ■ lO^^fc^A;^ + 
1.634 • lO^^fc^A;^ + 6.618 ■ lO^^kfkl - 4.437 ■ lO^^lk^ + 1.643 ■ lO^^t^t, 

representing the parameters where J(^4) has a pair of purely imaginary 
eigenvalues A3^4 = ±iujQ with 

0.0071(-2.6479 ■ 10^^ + A;2)(1.4219 ■ 10"^ + k2) 



ujQ = 1.2909 



0.6909 + 



A;2(3.1305- 10-7 + A;iA;2 
8.5745 ■ 10-7 / 1 



+ ..o... , uu.uu (1-0783 ■ 10- + (33) 



k2 V(6.2611- 10-7 + A;ifc2)fciA;2 

5.0825 ■ 10^%)A;2 - 3.1286 ■ 10-^(1.3649 ■ 10-^A;2) (6.0127 ■ 10"*^ + k2) + 



0.9391 A;?(9.5980 ■ 10"^ + A;2)(8.1563 ■ 10"*^ + fca 

13 



1/2- 



1/2 















/■J 

















> 



Figure 1: Set of admissible parameters 5 and Hopf curve S. 
Thus one has (see Fig. [1]) 

5 = 5+ US U5_. 

For parameter values in the region 5+ the equihbrium is unstable since 
the Jacobian matrix J{A4) has two complex eigenvalues with positive real 
parts and two other real negative eigenvalues. For parameter values in the 
region iS_ the equilibrium is asymptotically stable since J^Ai) has two 
complex eigenvalues with negative real parts and two other real negative 
eigenvalues. The curve S is the curve of admissible parameters where the 
equilibrium A4, is a Hopf point. 

Theorem 3.1 Consider the differential equations (QP with the parameters 
given in the table (E^. // (^1,^2) £ S then the two parameter family of 
differential equations (Op has a transversal Hopf point at A4. This Hopf 
point at A4 is unstable and for each {ki, G iS_, but close to S, there exists 
an unstable periodic orbit near the asymptotically stable equilibrium point A/^. 
See Fig\^ 
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Computer Assisted Proof. The proof follows the steps outlined in Sub- 
section I3.1[ However, all the expressions in the proof are too long to be put 
in print. For this reason, in the site [10] have been posted the main steps 
of the long calculations involved in the proof. This has been done in the 
form of a notebook for MATHEMATICA 5 [11]. A sufficient condition for 
being a Hopf point is that the first Lyapunov coefficient h ^ 0. In fact, it 
can be shown numerically that h{ki,k2) > for all values (fci,/c2) € S. A 
particular case and other related calculations are considered below for the 
sake of illustration. 

Take the particular point Q = {ki = 0.00331, fcs = 0.00100) G S with five 
decimal round-off coordinates [7J. For these values of the parameters 

Ai = (18543.57758,21634.17385,738.0525862,553.5394397). 

The Jacobian matrix J^A^) has eigenvalues 

Ai = -3.61058, Aa = -0.22912, A3,4 = ±2.84670i. 



/ 820.5542609 + 1080.774610Z \ 
295.1756045 - 139.5588184i 
^ ~ 862.8021803 + 130.4940530? ' 
\ 26.01486634 -87.2710071 7i / 

/ 0.00003314748646 + 0.00006274424412i \ 
-0.00003846764141 + 0.00003199241887? 
0.0005233172006 + 0.00007678211168i 
\ 0.001254520214 - 0.004888597529i / 



and thus 



ujo = 2.84670. 



(34) 



From ([T7|l the eigenvectors q and p have the form 



One has 



B{q,q) 



I -767.7418261 + 289.3499796i \ 


786.4902302 + 276.2661945? 
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and 



B{q,q) 



\ 



482.6477605 


-809.3875158 + 0.6 ■ 10" 




From (fT9|) and (120|) the complex vectors hu and /i2o have the form 



/ -1622.977370 + 0.9904140546- 10-^i \ 



-h 
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-1893.473598 + 0.1155483063 • 10 
-5.359116331 - 0.3117318364 ■ IQ-^ 
\ -4.019337253 - 0.2337988771 ■ IQ-^ ) 



/ 71.87520338 + 68.12253398i \ 
9.204672581 - 7.866965075i 
83.57142169 - 174.4554220i 
\ -8.839482676 - 5.024621730i / 

From (122|) the complex number G21 is given by 



h 



20 



G 



21 



0.057297 - 0.027485Z, 



(35) 



and from ( !23|) . ( !34l) and ( !35l) one has the first Lyapunov coefficient at Q 

h{Q) = 0.00353522 > 0. 



(36) 



The above calculations have also been checked with 10 decimals round-off 
precision performed using the software MATHEMATICA 5 [11]. See |10] . 

Some other values of pairs (^1,^2) € S, the values of the complex eigen- 
values of 7(^4) as well as the corresponding values of /i(/ci, ^2) are listed the 
table below. The calculations leading to these values can be found in [TO] . 
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hi 




r'OTTi'nlp'v pnxPTiA/'P liipc; of J{ JXa\ 




n nnn48i 

W . WW Wt:(J X (J 


000481 2 

W . W WWtiC ± ^ 


+4 764 ^ifi? 


4 69457 ■ 1 0~® 


W.WWW 1 UOrt 


0003535 

W.WWWOtJfJtJ 


+3 98051? 


1 21 597 ■ 1 0"'^ 


0.0011096 


0.0003086 


±3.5905H 


2.17614-10-^ 


0.0014238 


0.0002950 


±3.35518i 


3.29018 ■ 10-^ 


0.0017379 


0.0003001 


±3.19780i 


4.52683 ■ 10-^ 


0.0020521 


0.0003220 


±3.08560i 


5.86835 ■ 10-^ 


0.0023663 


0.0003649 


±3.00207i 


7.30346 ■ 10"'^ 


0.0026804 


0.0004427 


±2.93795i 


8.82421 ■ 10"^ 


0.0029946 


0.0005957 


±2.88762i 


1.04245 ■ 10-6 


0.0033088 


0.0009855 


±2.84745i 


1.20994-10-6 


0.0036230 


0.0035924 


±2.81501i 


1.38449 ■ 10-6 



Remark 3.2 The value of the first Lyapunov coefficient li does depend on 
the normalization of the eigenvectors q and p, while its sign is invariant 
under scaling of q and p obeying the relative normalization. See 19], p. 99. 
The values li{Q) in ^SE^ and li{ki,k2) in the above Table are obtained with 
two different choices of the eigenvectors q and p, see 110]. This explains the 
difference in the order of magnitude of the numbers involved. 

As a consequence of Theorem 3.1 there are no Hopf points of codimension 
2 on S since the sign of the first Lyapunov coefficient does not change. In 
Fig. [2] is illustrated the bifurcation diagram for a typical point on the curve 
S. 

Assuming the same values in the table fl29l) in next theorem we study 
the behavior of the Hopf curve S in the set of admissible parameters S (see 
equation (1321) ) as the parameter C2 increases. In fact, the carrying capacity, 
representing several other factors, has a determinant role on the populations 
under study. 

Theorem 3.3 The one parameter family of curves = A-^(O) has only 
one point of tangency T with the line ki = k2 for C2 = 650.41463. For values 
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C2 > 650.41463 the curve does not intersect the set S. Therefore for 
values 02 > 650.41463 the set S+ is empty, 5 = 5_ and the equilibrium A4 
is asymptotically stable for all values (fci, ^2) ^ S. See Fig. O 

Proof. The surface of Hopf points, or equivalently the one parameter family 
of Hopf curves, where J(^4) has a pair of purely imaginary eigenvalues is 
defined by = {A(A;i, k2, C2) = 0} where A{ki, k2, C2) (see ffTTl) ) is given by 

1699.422 - 2337.623c2A;i - 6.860 ■ 10^A;2 - 1114.941c^A;^ - 2.175 ■ lO^^A;^ - 
4.994 ■ 10^02^1^2 - 214.747c^A;? - 4.079 ■ 10^c^A;^A;2 - 1.529 • 10^^C2A;iA;^ - 
7.809 ■ 10^c^A;^A;2 - 4.319 ■ lO^^clklkl - 1.540 ■ 10^^C2A;iA;^ - 
4.755 • I0^c'^k% - 1.752 ■ lO^^c^fc^A;^ - 6.741 ■ lO^^c^A^^A;^ + 
2.940 ■ lO^^c^^ktkl - 1-703 ■ lO^^clkfkl + 1.634 ■ lO^^c^A^^A;^ + 
6.618 • lO^^ctktkl - 4.437 • lO^^clklk^ + 1.643 • lO^^ctk^kl 
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Figure 3: Curve S intersects S at one point T. 

The intersection of the surface with the plane ki = k2 determines the 
curve C, given imphcitly by 

N{ki, C2) = 1699.422 - (2337.623c2 + 6.860 ■ W)ki - (1114.941c^ + 
2.175 ■ 10^2 + 4.994 ■ lO^Ca)^;^ - (214.747c^ + 4.079 ■ lO^cj + 
1.529 ■ 10^^C2)kf - (7.809 ■ lO^c^ + 4.319 ■ lO^^c^ + 1.540 ■ 10^^C2)A;f - 
(4.755 • lO^^c^ + 1.752 • lO^^c^ + 6.741 ■ 10^^c^)A;i^ + (2.940 ■ 10^°c^ - 
1.703 ■ lO^^c^ + 1.634 ■ 10^2c^)A;^ + (6.618 ■ lO^^^c^ - 4.437 ■ \^^'^S^k\ + 

1.643 ■ lO^^c^A;^ = 0. 

Differentiating imphcitly the above expression with respect to k\ one has 

— £ = = 0, — J < 0, 

ac2 ^ 

at ki = 0.00035 and C2 = 650.41463. Therefore the curve C is a graph near 
the point {ki = 0.00035, C2 = 650.41463) and has a local maximum point at 
ki = 0.00035. It can be shown that this maximum is global since dc2/dki 
has no other zero. It is easy to verify through a calculation that the point 
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T = {ki,k2) = (0.00035,0.00035) belongs to S,^ for ca = 650.41463. Now 
the gradient of Ac^ at T for C2 = 650.41463 is given by 

(-1.73746 ■ 10^°, 1.73746 ■ 10^°), 

which is parallel to the vector (—1, 1), the normal to the line ki = ^2- 

■ 

Remark 3.4 Since does depend on the parameter C2, according to Eq. 

Ij^, so does the admissible region S = Sc2- For C2 = 650.41463, a calculation 
gives ki^^^ = 0.00338491. This is compatible with position ofT at ki = k2 = 
0.00035, as illustrated in Figure\^ 

4 Concluding Comments 

In this paper we studied the system ([1]) of interest as a mathematical model 
for biological control, proposed by Yang and Ternes [H El [2] and studied also 
by Santos [7]. Valuable field data are provided in [3], valid for the citrus 
leafminer and its native and imported enemies in the region of Limeira (Sao 
Paulo, Brazil). An extensive, enlightening discussion of the economic and 
agricultural interest of the problem, other pertinent differential equations 
models as well as extensive bibliography, are also presented there. 

Under conditions made explicit in Remark 12.11 we determine the unique 
equilibrium point (^4) with positive coordinates and establish necessary and 
sufficient conditions for its (Lyapunov) stability (Theorem 12.41) . It can be 
seen however that this condition A > 0, when expressed in terms of the 
parameters is a rational function whose denominator does not vanish and its 
numerator is a polynomial of too many terms to be put in print, but still 
amenable to numerical calculations. For this reason the treatment of the 
stability of (^4) in Subsection 13.21 is computer assisted. The conclusion of 
this study, made precise in Theorem 13.11 is the existence of periodic orbits 
obtained by Hopf bifurcation, on the side (of A = 0) where Aa is an attractor. 
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The study of the general analytic and geometric properties of the bound- 
ary of the stability region, given by the Hopf variety A = 0, so as to include 
parameter values of biological interest as proposed here as well as others 
appearing in the work of Ternes [3], remain at the present moment as a 
mathematical challenge. Theorem 13.31 gives only a thin slice of the geometry. 

The reports in [12] and [13], among many others, show that the interest 
for the combat of the citrus leafminer extends to most regions where citrus 
trees grow. 

The Mathematica notebooks [10], with the table (129|) . used in the com- 
puter assisted arguments for the proofs of Theorems 13.11 and 13. 3[ can be 
adapted to tables with data pertinent to other geographic and climatic re- 
gions and involving different host-parasitoid interactions. 

In [2] Ternes and Yang discuss, with pertinent documentation, the intro- 
duction of a foreign parasitoid, Ageniaspis citricola to add the native Galeop- 
somya Fausta in the combat with the leafminer, Phyllocnistis citrella. They 
propose a model with eight differential equations for the three species and 
their immature stages. In [2] and [3] are given starting steps for an analysis 
of the stability of the equilibria in this extended eight - dimensional system. 
Based in a numerical study of a complex equilibrium point they recommend 
that the biological control of the leafminer he implemented with both the na- 
tive and foreign parasitoids. Meanwhile, the Hopf bifurcation analytic and 
computer algebra study of the complex equilibria of the eight equations, with 
the methods used in the present paper, seems unsurmountable at the present 
moment, due to the large number of parameters involved. 
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